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Abstract 

In this paper, we propose simple estimation methods dedicated to a semiparametric family 
of bivariate copulas. These copulas can be simply estimated through the estimation of their 
univariate generating function. We take proht of this result to estimate the associated mea- 
sures of association as well as the high probability regions of the copula. These procedures 
are illustrated on simulations and on real data. 

Keywords: Copulas, nonparametric estimate, measures of association, high probability 
regions. 
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1 Introduction 



The theory of copulas provides a relevant tool to build multivariate probability laws, from fixed 
margins and required degree of dependence. From Sklar's Theorem [22], the dependence prop- 
erties of a continuous multivariate distribution H can be entirely summarized, independently 
of its margins, by a copula, uniquely associated with H. Several families of copulas, such as 
Archimedian copulas [8] or copulas with polynomial sections [2Q|, I18j have been proposed. More 
recently, we proposed to give up the polynomial form to work with a semiparametric family of 
copulas pp. This permits to increase the dependence degree and to preserve the dependence 
properties of copulas with polynomial sections without significantly complexifying the model. 
Furthermore, the family of copulas is generated as simply as Archimedian copulas, that is by an 
univariate function. 

Among the numerous papers dedicated to the construction of copulas, very few of them pro- 
pose some associated inference procedures. Besides, starting from data, it is very difficult to 
find "Which copula is the right one?", see |12j for a review on this problem in the financial 
modeling context. The first attempt to estimate copulas is achieved in [2] with the introduction 
of nonparametric estimates based on empirical copulas. A recent contribution to the nonpara- 
metric estimation of copula for time series is presented in [23]. An alternative approach is to 
choose a parametric family of copulas and to estimate the parameters by a maximum likelihood 
method [15]. Both approaches suffer from their own drawbacks. In the first case, a fully non- 
parametric estimate is likely to suffer from the so-called "curse of dimensionality" . It would have 
a high variance for large number of margins and moderate size of datasets. At the opposite, a 
parametric estimate can lead to a very high bias if the prior family of copulas is not appropriate. 
Restricting to Archimedian copulas, Genest &: Rivest [9] proposed a semiparametric estimate. In 
this family, estimating the copula reduces to estimating the generating function. It is therefore 
realistic to consider a nonparametric estimate of this univariate function. Here, we show that, 
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similarly, the estimation of a copula in the semiparametric family p] can be simply achieved by 
estimating nonparametrically the univariate generating function. This permits to overcome the 
curse of dimensionality. We deduce of this estimator some statistical procedures for estimating 
the associated dependence coefficients and high probability regions. 

In section [21 the expression of the considered semiparametric family of copulas is given and its 
basic properties are recalled. Section [3] is dedicated to the estimation of the generating function 
and of the dependence coefficients. In section [H we address the problem of estimating high prob- 
ability regions. All the proposed estimates are experimented on simulated samples in section [5] 
and on real data in section [U 

2 Definition and basic properties 

Throughout this paper, we note I = [0, 1]. A bivariate copula defined on the unit square I 2 is 
a bivariate cumulative distribution function with univariate uniform I margins. Equivalently, it 
must satisfy the following properties : 

(PI) C{u,0) = C(0,v) = 0, V(u,v) G I 2 , 

(P2) C{u, l) = u and C(l, v) = v, V(u, v) G I 2 , 

(P3) A(m, u 2 , vt, v 2 ) = C(u 2 ,v 2 ) - C(u 2 , vt) - C(u\,v 2 ) + C(ui,vi) > 0, V(ui, u 2 , v 1 ,v 2 ) £ I 4 , 
such that u\ < u 2 and v\ < v 2 . 

Let us recall that, from Sklar's Theorem [22], any bivariate distribution with cumulative distri- 
bution function H(x, y) and marginal cumulative distribution functions F{x) and G(y) can be 
written H(x,y) = C(F(x),G(y)), where C is a copula. This result justifies the use of copulas 
for building bivariate distributions. 

Here, we consider the semiparametric family of functions defined on I 2 by: 

Ce^u, v)=uv + 6<f>(u)<f>(v), 6d [—1,1], (2.1) 
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where <f> is a function on I. This family was first introduced in [21], chapter 3, and is a particular 
case of Farlie's family [6j. It is extensively studied in pQ. In particular, the following basic lemma 
is proved: 

Lemma 1 <j) generates a parametric family of copulas {Cq^, 6 G [—1,1]} if and only if it 
satisfies the following conditions : 



The function <f> plays a role similar to the generating function in Archimedian copulas [8J . Each 
copula Cq a is entirely described by the univariate function <f> and the parameter 6, which tunes 
the dependence between the margins. For instance, cp(x) = x(l—x) generates the Farlie-Gumbel- 
Morgenstern (FGM) family of copulas [T7], which contains all copulas with both horizontal and 
vertical quadratic sections |20| . Another example is (j>(x) = x(l — x)(l — 2x) which defines the 
parametric family of symmetric copulas with cubic sections proposed in [18], equation (4.4). Of 
course, it is also possible to choose <p so as to define new copulas, see section [5] for an example. 

3 Estimation of the generating function 
3.1 Preliminaries 

Let (X, Y) be a random pair from the cumulative distribution function H{x, y) = Cq < k(F(x), G(y)), 
where F(x) and G(y) are respectively the cumulative distribution functions of X and Y. The 
estimation of the copula reduces to estimating the generating function <ft and the parameter 9. 
This estimation clearly suffers from an identifiability problem since, for instance, replacing <f> 
by 4>/\fa and 6 by a.0 for any positive a leads to the same copula. When 0^0, introducing 
s = 0/\6\ and ip = \/\6\(j) yields 




(i) 0(0) = <P(1) = 



(ii) \4>(x) - 4>(y)\ < \x - y\ for all (x,y) G I 2 . 




Ce,<f>(u, v) = C s ^(u, v) = uv + sip(u)ip(v) 



(3.1) 
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where ip satisfies the conditions of Lemma [T] and s G { — 1,1}. The identifiability problem is 
not fully overcome with the new parameterization (13. ip since the sign of ip cannot be identified. 
Thus, we limit ourselves to the Positively Quadrant Dependent (PDQ) context. Recall that X 
and Y are PQD [14] . section 2.1.1, if and only if 

V(s,y) G M 2 , P(X < x,Y < y) > P(X < x)P(Y < y). 
As shown in pQ, theorem 3, X and Y are PQD if and only if 

6 > and, either Vu G J, 0(it) > or Vu G /, 4>{u) < 0. (3.2) 
If (X,Y) are PQD, then the copula (|3.ip can always be rewritten as: 

Citf(u, v) = uv + ip(u)ip(v), (3.3) 

where ^ is a non negative function satisfying the conditions of Lemma [TJ In the following, we 
limit ourselves to the estimation of C±^, or equivalently to the estimation of ip in this context. 
We refer to section [7] for possible improvements of the estimation method. 

3.2 Estimation of ip(w) 

Let {(x{, yi), i = 1, . . . , n} a sample of (X, Y) from the cumulative distribution function H(x, y). 
The rank transformations U{ = Rank(xj)/n and Vi = Rank(yj)/n yield an approximate sample 
from the copula C\^{u,v). The estimation of ip(w) relies on the pseudo-observations Wi = 
max(uj, Vi), i = 1, . . . , n which have the common distribution function C\^(w, w) = w 2 + ip 2 (w). 
In order to obtain a regular estimated function, this estimate is written as a linear combination 
of a denombrable set A of functions: 

#"•0 = a k£k{w). (3.4) 
keA 

The set of functions {ek{w), w £ I, k G A} need not to be orthogonal but the condition 
efc(O) = efc(l) = is required for all k G A in order to ensure ip(0) = ip(l) = 0, and thus to 
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respect condition (i) of LemmaCQ Introducing w\ jTl < ■ ■ ■ < w njTl the ordered statistics associated 
to the Wi, i = 1, . . . , n, the coefficients ak are determined by the constrained least-square problem 

a = argmin{|| Ma - b\\ 2 , < (Ma),, -1 < (M'a)i < 1, i = l,...,n} (3.5) 

where M and M' are two matrices such that M^k = ek(wi in ) , M' ik = e' k {wi^ n ) for k G A, 
i G {1, . . . ,n}, and b is a vector defined by bi = (i/(n + 1) — ^ 2 n ) 1//2 , i G {1, ... The 
minimization of || Ma — 6|| 2 ensures that for all i = 1, . . . , n 

■4>{wi,n) 2 = Ci t ^(w i>n , w iyn ) - wf n wi/(rt + l)- wf n . (3.6) 

The positivity conditions < (Ma)i impose that < ip(wi tn ), and the bound conditions 

— 1 < (M'a)i < 1 are interpreted as (wi jn )\ < 1 which allows to fulfil Lemma QJii). 

In practice, the optimization process provides reasonably sparse solutions, that is vectors a with 

a moderate number of nonzero components. 

3.3 Estimation of an association measure 

Two measures of association between the components of the random pair (X, Y) are usually 
considered |14| . section 2.1.9. The Kendall's Tau is the probability of concordance minus the 
probability of discordance of two random pairs {X\,Y\) and (X2, Y2) described by the same joint 
bivariate law H(x,y) = C{F(x),G{y)). It only depends on the copula: 

r = 4/ C(u,v)dC(u,v)-l. (3.7) 
Jp 

The Spearman's Rho is the probability of concordance minus the probability of discordance of 
two pairs {X\,Y-\) and (^2,12) with respective joint cumulative law H(x,y) and F(x)G(y), 

p = 12 C(u,v)dudv — 3. 
Jp 

In PP, proposition 1, it is shown that, within the family (|2.ip . these two measures are equivalent 
up to a scale factor: 2p = 3r. Here, we focus on the Spearman's Rho whose expression is in our 
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family: 

p = 126 (J 4>(u)du\ = 12 (J i>(u)du 
We then propose an estimate based on the estimate ip: 

p SP = 12(j2 a ^k] > ( 3 - 8 ) 
\keA J 

where we have introduced (5k = fr e k( u )du- This integral can be calculated analytically or 

numerically by Simpson's rule, depending on the complexity of the basis of functions. It is 

also possible to estimate p in a nonparametric way by rescaling the empirical version of (|3.7p 

introduced in [9] with the factor 3/2 to obtain: 

6 n n 3 
= n ( n -l) S X K < u » v i < v d ~ o ' 

where 1{.} is the indicator function. The two estimates p SP and p NP are compared on simulations 
in section [5j 

4 Estimation of high probability regions 
4.1 The general problem 

Let us recall the definition of a p-dimensional a-quantile of a distribution P. Let S the class of 
Borel measurable sets of W and let A be the Lebesgue measure defined on S: 

Q a = inf{A(5) : P{S) > a, S € S}, < a < 1. 

Here, Q a is the minimum volume S E S that contains at least a fraction a of the probability 
mass. This is a particular case of the general quantile function introduced by Einmal and 
Mason [5j. A particular attention has been paid to the estimation of Q\, the support of the 
distribution, from a sample {Mi, . . . ,M n } of M p . The early paper of Geffroy [7] takes place in 
the case p = 2 and the considered supports are written 

Ql = {(x, y)eM. 2 :0<x<l; 0<y< f(x)}, 
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where / is an unknown function. More recently, smooth estimates of the frontier function / have 
been proposed [H] and their extension to star-shaped supports has been studied [13]. Numerous 
works have been dedicated to the nonparametric estimate 

n 

Ql = {jB(Mi,r n ), 
i=i 

where B(Mi,r n ) is the ball of radius r n and centered at Mj, see for instance [3] [TO]. In the latter 
paper, a partition based estimate is also proposed. Introducing {Ki}i>± a partition of R 2 , the 
estimate is simply: 

Q 1 = }jK i l{M j £Ki}. 

In subsection 14.21 we propose an estimate of Q a , a < 1, when P is a bivariate copula, with 
a similar principle. Subsection 14.31 examines the special case of the semiparametric family of 
copulas {Ci^}. The estimation of such high probability regions is of particular interest for 
bivariate copulas since some dependence properties can be read on this graph, see [16]. To this 
end, introduce the two diagonal lines £\\ v = u and £2- v = 1 — u of the unit square I 2 . If 
the graph of Q a is nearly symmetric with respect to both diagonals l\ and £2 then the copula 
models weakly dependent variates. On the contrary, if the graph is concentrated along l\, then 
the copula models strongly positively dependent variates. 

4.2 The algorithm 

Let {Ifc, k = 1, . . . , N} be the equidistant iV-partition of / and K^^ = 1^x1^ the associated 
N x N grid. Denote 6k,e £ {0, 1}, k = 1, . . . , N, £ = 1, . . . , N a binary N x N array and introduce 
the estimate 

Q a = \jK k/ l{S k ,i = l}, (4.1) 

k,e 

where the 5k,e are defined by the optimization problem 

. JV N 

k=l 1=1 
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under the constraints 6 k g S {0, 1} and 

N N 

^^<5 M P(if M )>a. (4.3) 

k=l 1=1 

The quantity P(K k g) is an estimation of the probability P{K k g). The quality of the estimate 
Q a strongly depends on the quality of this estimate. This is discussed in subsection 14.31 Let 
us note that (|4.2[) is equivalent to minimizing X(Q a ) and (j4.3|) corresponds to the constraint 
P{Qa) > a - This optimization problem can be solved with a simple algorithm. The first step 
consists to sorting the estimated probabilities P(K k e) in decreasing order to obtain the sequence 
P T , t = 1, . . . , TV 2 . The second step is the computation of the number of subsets of the partition 
which are going to be used: 

J = min jj, y~]Pr > « j • 
The last step is the selection of the J first subsets: 8 k t = 1 if 1 < r(k,£) < J, which leads to 
the estimate (j4.1j) . 

4.3 Estimation of P(K k) e) 

If no information is available on the distribution P, one can use the nonparametric estimate 

1 - 

PMK k ,e) = -JZHMi €K k4 }. 
n 

i=l 

Now, restricting ourselves to the family {Ci,/,} allows to consider semiparametric estimates, 
which are much more accurate than the nonparametric one. Taking into account of (|3.3f) . it 
follows that 

t>( V ^ ^ ( k l \ „ ( k ~ l *\ n ( k n f k ~ l 

P{K k/ ) - j+dj^— ,— 

1 U^^-^^Ud)-^"- 1 



N 2 \ \Nj V N J J \ \ N J \ N 
Thus, basing on section El the following semiparametric estimate can be introduced 

The results using P SP and P NP are compared on simulations in the next section. 
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5 Simulation results 

All the numerical experiments of this section have been conducted on the family of copulas 
generated by the set of functions 

Vfc>l, ip k (x) = l- (x fc + (l-x) fc ) lA , xel. (5.1) 

For the sake of simplicity, the resulting family of copulas will be denoted by Ck(u, v) = Cx^Aii, v). 
This family is interesting since it can represent very different distributions: 

• When k = 1, Ci(x,y) is the uniform distribution on the unit square I 2 . The associated 
Spearman's Rho is p\ = 0. 

• Letting k — > oo, we obtain ipk(x) — > ipoo(x) = min(x, 1 — x) for all x £ I. Thus, it 
appears that C 00 (x,y) is a mixture of two uniform distributions on the squares [0, 1/2] 2 
and [1/2, l] 2 with mixing parameter 1/2. The associated Spearman's Rho is p^ = 3/4, 
the maximum value in the family (12. ip . 

• When 1 < k < oo, we get a bivariate distribution "interpolating" between the two previous 
one. Up to our knowledge, it is not possible to calculate explicitly. 

5.1 Simulation of data from the copula 

The simulation method described in [19], p. 36 is used. First, two independent uniform samples 
Ui and ti, i = 1, . . . , n are simulated. Second, for each % = 1, . . . , n, fj is computed such that 

d ^ , . 
ou 

by a dichotomy procedure. The resulting sample (ui,Vi), i = l,...,n has joint distribution 
function Ck(u,v). 
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5.2 Estimation of the generating function 

In this paragraph, we have chosen n = 100. Starting from the sample (ui,Vi), % = 1, ... ,n, the 
generating function is estimated using the procedure described in section [3j The chosen basis 
of functions is doubly-indexed by a scale parameter s and a location parameter I: 



where A = {(s,£), < i < 2(2 S — 1), < s}. See figure[I]for a graph of the first basis functions. 
In figure[2j the estimation of ip2(x), t/JA (x) and ips( x ) is compared to the true generating functions. 
These first results are visually satisfying. The optimization procedure selects about 30 basis 
functions. A more precise comparison is achieved in table Q] by repeating each estimation on 
100 different samples. On the basis of these 100 repetitions, the mean value and the standard 
deviation of the L2 error 



are evaluated, as well as the mean value and the standard deviation of the two estimates p SP 
and p NP of the Spearman's rho. The integral appearing in p^ is calculated numerically using 
Simpson's rule. The computation of p SP is based on (I3.8P and is thus explicit after remarking 
that S> £ = fj e Sj i(x)dx = 2 l ~ s /ir. 

The mean values of e (about 10 -2 ) confirm that the function ip^ is correctly estimated. Moreover 
it appears that the semi parametric estimation of the Spearman's Rho is better than the non 
parametric estimate, excepted for the case k = 1. The standard deviations are similar for the 
two estimates. 

5.3 Estimation of high probability regions 

Starting from the estimations of ip2(x), ipi(x) and ips(x) obtained above, it is possible to es- 
timate high probability regions Q a with the procedure described in section The following 
probabilities a = 0.25, a = 0.5 and a = 0.75 are considered. Here, the regions obtained using 




(5.2) 




1/2 



£ = 
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the semiparametric estimate P SP , the nonparametric estimate P NP and the true probability P 
can be compared. Of course, this probability depends on ip(x) and cannot be used in practical 
situations. Here n = 500 and the estimated regions are obtained by a discretisation on a N x N 
grid with N = 30 when using P SP or P and N = 8 when using P NP so as to obtain approxima- 
tively 500/64 ~ 8 points in each cell. The results are presented on figures [3H5J It is apparent 
that the estimations obtained with the true probability P and its semiparametric estimate P SP 
are very close, especially for moderate values of k. This confirms the results obtained in the 
previous paragraph. On the contrary, the nonparametric estimate accuracy is very poor for such 
values of the sample size n. We can also observe that, as k increases, the high probability regions 
are more and more concentrated in the neighborhood of the l\ diagonal line. This, together with 
table [H illustrates the property that the positive dependence is increasing with k. 



6 Real data 

The dataset consists of n = 225 countries on which two variables have been measured: X, the 
life expectancy at birth (years) in 2002 of the total population and Y, the difference between 
the life expectancy at birth of women and men. The data is available on the following web-site: 



http : //www. ode i .gov/ cia/publications/f actbook/ According to the PQD test proposed in 



these data are PQD. The first step is to compute the rank transformations Ui = Rank(xj)/n 
and Vi = Rank(yj)/n, see figure [7J Then, the generating function ip(x) is estimated using 
the basis (|5.2p . The optimization procedure yields an expansion of ip(x) with respect to 24 
basis functions. The function ip(x) is plotted in figure [6l The estimated Spearman's rho are 
p NP = 52.4% and p SP = 40.7%. These values seem to confirm a moderate positive dependence 
between the two variables. The higher the life expectancy is, the more important the difference 
between women and men is. Finally, one can estimate the high probability regions. The same 
parameters as in the simulation section are used. The results are presented in figure[7J Of course, 
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since the true function ip(x) is not known (and perhaps does not exist), it is only possible to 
compare the estimations obtained with the nonparametric and semiparametric estimates P NP 
and P SP . The semiparametric estimation reveals two main groups of countries. In the first 
one, both men and women share a small life expectancy with weak differences between sexes. 
In the second one, both men and women share a large life expectancy but with important 
differences between sexes. The sample size is not large enough to obtain meaningful results with 
the nonparametric estimate. 

7 Further work 

We have presented a method for estimating copulas in the bivariate family of copulas C s ^(u, v) 
in the PQD case. Even though a test [3j [23] has not rejected the PQD assumption, deciding if 
the copula model C s ^{u,v) is adapted to a particular dataset is an opened problem. It could 
be possible to build a goodness-of-fit test based on the comparison of the estimations p NP and 
p SP . The test would reject the model C s ^(u,v) if the difference \p NP — p SP \ is too large. If the 
PQD assumption is rejected, the proposed estimation method cannot be used. To overcome the 
identifiability problem, a possible modification of the method would be to replace the projection 
step (|3.4j) by the selection of the function tp in a database leading to the best approximation (I3.6p . 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 1: Graph of the first basis functions of (|5.2I) , Solid line: eofi(x), dotted line: ei t o(x), 
dashdot line: ei i i(x), dashed line: e\^{x). 



k 


Pk 


mean( / 5 SP ) 


std(/5 SP ) 


mean(/5 NP ) 


std(/3 NP ) 


mean(e) 


std(e) 




xlO" 2 


xlO" 2 


xl0~ 2 


xl0~ 2 


xlO" 2 


xlO" 2 


xl0~ 2 


1 





0.81 


(6.62) 


0.18 


(11.0) 


8.75 


(3.58) 


2 


42.5 


43.0 


(9.91) 


41.2 


(9.H) 


3.67 


(1.26) 


4 


66.4 


65.8 


(7.11) 


64.3 


(6.11) 


3.01 


(1.33)) 


6 


71.2 


70.6 


(7.60) 


68.8 


(6.09) 


3.10 


(1.35) 


8 


72.8 


72.1 


(7.68) 


70.2 


(6.05) 


3.10 


(1.17) 



Table 1: Estimation of the generating function and of the Spearman's Rho (pk)- The mean 
value and the standard deviation of the L2 error e as well as of the estimates p SP and p NP are 
evaluated on 100 repetitions. 
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(a) True generating functions ipk( x ), k G {2,4,8} 



0.5 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(b) Estimated generating functions ipk(x), k E {2,4,8} 



Figure 2: Comparison of the true generating functions ipk( x ) denned by (|5.ip and the estimated 
ones ipk(x). Solid line: k = 2, dashed line: k = 4, dashdot line: k = 8. 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 




Figure 3: Estimation of high probability regions Q a from C^m, v). Dark Grey: a = 0.25, grey: 
a = 0.5, light grey: a = 0.75. Top left: simulated sample, top right: nonparametric estimate, 
bottom left: semiparametric estimate, bottom right: semiparametric estimate with the true 
function ip. 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 



Figure 4: Estimation of high probability regions Q a from C&(u, v). Dark Grey: a = 0.25, grey: 
a = 0.5, light grey: a = 0.75. Top left: simulated sample, top right: nonparametric estimate, 
bottom left: semiparametric estimate, bottom right: semiparametric estimate with the true 
function ip. 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 



Figure 5: Estimation of high probability regions Q a from Cs(u,v). Dark Grey: a = 0.25, grey: 
a = 0.5, light grey: a = 0.75. Top left: simulated sample, top right: nonparametric estimate, 
bottom left: semiparametric estimate, bottom right: semiparametric estimate with the true 
function ip. 
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Figure 7: Estimation of high probability regions Q a from real data. Dark Grey: a = 0.25, grey: 
a = 0.5, light grey: a = 0.75. Top left: real data, top right: real data after rank transformation, 
bottom left: nonparametric estimate, bottom right: semiparametric estimate. 
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